Studies 1-4: Figures

Luhrmann, Weisman, et al.

In this notebook, we create figures for Studies 1-4.

source("./scripts_general/dependencies.R")

Loading required package: piecewiseSEM
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'lme4':
  method                          from
  cooks.distance.influence.merMod car 
  influence.merMod                car 
  dfbeta.influence.merMod         car 
  dfbetas.influence.merMod        car 

  This is piecewiseSEM version 2.0.2

  If you have used the package before, it is strongly recommended you read Section 3 of the vignette('piecewiseSEM') to familiarize yourself with the new syntax

  Questions or bugs can be addressed to <jlefcheck@bigelow.org>
Loading required package: lme4
Loading required package: Matrix
Loading required package: lmerTest
Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

Loading required package: brms
Loading required package: Rcpp
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'xts':
  method     from
  as.zoo.xts zoo 
Loading 'brms' package (version 2.10.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: ‘brms’

The following object is masked from ‘package:lme4’:

    ngrps

Loading required package: psych

Attaching package: ‘psych’

The following object is masked from ‘package:brms’:

    cs

Loading required package: kableExtra
Loading required package: langcog

Attaching package: ‘langcog’

The following object is masked from ‘package:base’:

    scale

Loading required package: tidyverse
── Attaching packages ─────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.2.1     ✔ purrr   0.3.3
✔ tibble  2.1.3     ✔ dplyr   0.8.3
✔ tidyr   1.0.0     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0
── Conflicts ────────────────────────────────────────────── tidyverse_conflicts() ──
✖ ggplot2::%+%()      masks psych::%+%()
✖ ggplot2::alpha()    masks psych::alpha()
✖ tidyr::expand()     masks Matrix::expand()
✖ dplyr::filter()     masks stats::filter()
✖ dplyr::group_rows() masks kableExtra::group_rows()
✖ dplyr::lag()        masks stats::lag()
✖ tidyr::pack()       masks Matrix::pack()
✖ tidyr::unpack()     masks Matrix::unpack()
Loading required package: cowplot

********************************************************
Note: As of version 1.0.0, cowplot does not change the
  default ggplot2 theme anymore. To recover the previous
  behavior, execute:
  theme_set(theme_cowplot())
********************************************************

Loading required package: sjstats

Attaching package: ‘sjstats’

The following objects are masked from ‘package:psych’:

    pca, phi

source("./scripts_general/custom_funs.R")
source("./scripts_general/var_recode_contrast.R")
source("./study1/scripts_s1/s1_var_groups.R")
source("./study2/scripts_s2/s2_var_groups.R")
source("./study3/scripts_s3/s3_var_groups.R")
source("./study4/scripts_s4/s4_var_groups.R")

setwd("./study1/analysis/")

The working directory was changed to /Users/anthrouser/Documents/Projects/Mind and Spirit overview paper/mindspiritquant/study1/analysis inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

source("../../scripts_general/data_load.R")

Parsed with column specification:
cols(
  .default = col_double(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  site = col_character(),
  religion = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character(),
  study = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  date = col_date(format = ""),
  researcher = col_character(),
  country = col_character(),
  quad = col_character(),
  subject_gender = col_character(),
  subject_job = col_character(),
  subject_schedule = col_character(),
  subject_livedhere = col_character(),
  subject_lang = col_character(),
  subject_marital = col_character(),
  subject_hs = col_character(),
  subject_school = col_character(),
  subject_liveswith = col_character(),
  servicesperweek = col_character(),
  godexpviaawe_freq = col_character(),
  sleephabit = col_character(),
  researcher_date = col_date(format = ""),
  researcher_notes = col_character(),
  site = col_character(),
  religion = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_double(),
  ctry = col_character(),
  wher = col_character(),
  recr = col_character(),
  whoc = col_character(),
  demo_chur = col_character(),
  demo_ethn = col_character(),
  demo_maj = col_character(),
  demo_pocc = col_character(),
  demo_rlgn = col_character(),
  demo_sex = col_character()
)
See spec(...) for full column specifications.
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_character(),
  X1 = col_double(),
  epi_subj = col_double(),
  epi_demo_age = col_double(),
  epi_demo_ses_num = col_double(),
  epi_demo_howr_num = col_double(),
  epi_demo_affr_num = col_double(),
  epi_demo_tung_num = col_double(),
  epi_demo_ytng_num = col_double(),
  epi_demo_affr_cat = col_logical(),
  epi_demo_tung_cat = col_logical(),
  epi_demo_ytng_cat = col_logical(),
  epi_version = col_double(),
  epi_charc = col_double()
)
See spec(...) for full column specifications.
2 parsing failures.
 row       col expected  actual                           file
1025 epi_charc a double Unclear '../../study3/data/d_demo.csv'
1027 epi_charc a double Unclear '../../study3/data/d_demo.csv'
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  score = col_double()
)
Joining, by = c("epi_ctry", "epi_subj")
funs() is soft deprecated as of dplyr 0.8.0
Please use a list of either functions or lambdas: 

  # Simple named list: 
  list(mean = mean, median = median)

  # Auto named with `tibble::lst()`: 
  tibble::lst(mean, median)

  # Using lambdas
  list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
This warning is displayed once per session.Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  epi_ctry = col_character(),
  epi_subj = col_double(),
  question = col_character(),
  response = col_double(),
  order = col_double(),
  question_text = col_character()
)
Joining, by = c("epi_ctry", "epi_subj", "question", "response", "order", "question_text")
Joining, by = c("epi_ctry", "epi_subj")
Column `epi_ctry` joining character vector and factor, coercing into character vectorJoining, by = "epi_subj"
Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_double(),
  p7_ctry = col_character(),
  p7_abs_check = col_character(),
  p7_dse_check = col_character(),
  p7_se_check = col_character(),
  p7_unev_check = col_character(),
  p7_exsen_check = col_character(),
  p7_por_check = col_character(),
  p7_mm_check = col_character(),
  p7_dem_sex = col_character(),
  p7_dem_pocc = col_character(),
  p7_dem_major = col_character(),
  p7_dem_ethnicity = col_character(),
  p7_dem_rur.urb = col_character(),
  p7_dem_affrd.basics = col_character(),
  p7_dem_religion = col_character(),
  p7_dem_church = col_character(),
  p7_dem_holy.tung.gif = col_character(),
  p7_abs_child.exp_cat = col_logical(),
  p7_abs_poetic_cat = col_logical(),
  p7_abs_tv.real_cat = col_logical()
  # ... with 162 more columns
)
See spec(...) for full column specifications.
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("p7_ctry", "p7_subj", "question", "response", "scale")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")
Joining, by = c("study", "p7_ctry", "p7_subj", "abs_score", "cog_score", "ctl_score", "dse_score", "hall_score", "para_score", "por_score", "pv_score", "spev_score")

Summary statistics

Spiritual Events scores

Study 1

d1 %>%
  ggplot(aes(x = country, y = spirit_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spirit_score, na.rm = T),
                              sd = sd(spirit_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Study 1: Spiritual Events score (range: 0-1)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 2

d2 %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 2: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 3

d3 %>%
  # correct for scaling in original dataset
  mutate(spirit_score = spirit_score * 4) %>%
  ggplot(aes(x = epi_ctry, y = spirit_score, color = epi_ctry, fill = epi_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(epi_ctry) %>%
                    summarise(mean = mean(spirit_score, na.rm = T),
                              sd = sd(spirit_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 3: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 4

d4 %>%
  ggplot(aes(x = p7_ctry, y = spev_score, color = p7_ctry, fill = p7_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(p7_ctry) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 4: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

All studies

d_all <- d1 %>%
  select(study, country, site, religion, subject_id, 
         por_score, abs_score, spirit_score) %>%
  rename(spev_score = spirit_score,
         pv_score = por_score) %>%
  mutate(religion = recode_factor(religion,
                                  "indigenous" = "Indigenous Religion",
                                  "charismatic" = "Charismatic Christianity"),
         # rescale to 0-1
         pv_score = pv_score/3) %>%
  full_join(d2 %>% 
              select(study, country, subj, 
                     abs_score, spev_score, dse_score) %>% 
              rename(subject_id = subj) %>%
              # rescale to 0-1
              mutate(spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  full_join(d3 %>% 
              select(study, epi_ctry, epi_sample, epi_subj, 
                     por_score, spirit_score) %>%
              rename(country = epi_ctry, 
                     religion = epi_sample,
                     subject_id = epi_subj,
                     spev_score = spirit_score) %>%
              mutate(religion = recode_factor(religion,
                                              "general population" = "General Population",
                                              "charismatic" = "Charismatic Christianity"))) %>%
  full_join(d4 %>%
              select(study, p7_ctry, p7_subj, 
                     por_score, pv_score, abs_score, spev_score, dse_score) %>%
              rename(country = p7_ctry,
                     subject_id = p7_subj) %>%
              # rescale to 0-1
              mutate(por_score = por_score/2,
                     pv_score = pv_score/3,
                     spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  mutate(religion = factor(religion, 
                           levels = c("General Population", 
                                      "Indigenous Religion",
                                      "Charismatic Christianity")),
         study = gsub("study", "Study", study)) 

Joining, by = c("study", "country", "religion", "subject_id", "abs_score", "spev_score")
Column `religion` joining factor and character vector, coercing into character vectorJoining, by = c("study", "country", "religion", "subject_id", "spev_score")
Column `religion` joining character vector and factor, coercing into character vectorJoining, by = c("study", "country", "religion", "subject_id", "pv_score", "abs_score", "spev_score", "dse_score", "por_score")

d_all %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country,
             group = study)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.25,
                                             jitter.height = 0,
                                             dodge.width = 0.75), 
             alpha = 0.1) +
  geom_pointrange(data = d_all %>%
                    group_by(study, country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, shape = study), 
                  position = position_dodge(width = 0.75),
                  color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom") +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       shape = "Study", 
       caption = "Error bars are ±1 standard deviation from the mean")

d_spev_sum <- d_all %>%
  filter(!is.na(spev_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(spev_score, na.rm = T),
            sd = sd(spev_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>%
  ggplot(aes(x = country, y = spev_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_spev_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_spev_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Porosity scores

d_por_sum <- d_all %>%
  gather(por_scale, por_score, c(pv_score, por_score)) %>%
  mutate(por_scale = recode(por_scale,
                            "pv_score" = "Porosity Vignettes",
                            "por_score" = "Porosity Scale")) %>%
  filter(!is.na(por_score)) %>%
  group_by(study, country, religion, por_scale) %>%
  summarise(mean = mean(por_score, na.rm = T),
            sd = sd(por_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>% 
  gather(por_scale, por_score, c(pv_score, por_score)) %>%
  mutate(por_scale = recode(por_scale,
                            "pv_score" = "Porosity Vignettes",
                            "por_score" = "Porosity Scale")) %>%
  filter(!is.na(por_score)) %>%
  ggplot(aes(x = country, y = por_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(cols = vars(study, por_scale), scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_por_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_por_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Porosity score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Absorption scores

d_abs_sum <- d_all %>%
  filter(!is.na(abs_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(abs_score, na.rm = T),
            sd = sd(abs_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>% 
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = country, y = abs_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(. ~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_abs_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_abs_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Absorption score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

r1_spev <- lm(spirit_score_std ~ country, d1)
r1_pv <- lm(por_score_std ~ country, d1)
r1_abs <- lm(abs_score_std ~ country, d1)

r2_spev <- lm(spev_score_std ~ country, d2)
r2_dse <- lm(dse_score_std ~ country, d2)
r2_abs <- lm(abs_score_std ~ country, d2)

r3_spev <- lm(spirit_score_std ~ epi_ctry, d3)
r3_por <- lm(por_score_std ~ epi_ctry, d3)

r4_spev <- lm(spev_score_std ~ p7_ctry, d4)
r4_dse <- lm(dse_score_std ~ p7_ctry, d4)
r4_dse <- lm(dse_score_std ~ p7_ctry, d4)
r4_por <- lm(por_score_std ~ p7_ctry, d4)
r4_pv <- lm(pv_score_std ~ p7_ctry, d4)
r4_abs <- lm(abs_score_std ~ p7_ctry, d4)

var	scale	study 1	study 2	study 3	study 4
spiritual experience	DSE	.	.	50%	50%
	spiritual events	22%	24%	38%	31%
porosity	porosity scale	.	61%	.	39%
	porosity vignettes	37%	.	.	27%
absorption	absorption	10%	.	11%	9%

r1_spev_pv <- lm(spirit_score_std ~ por_score_std, d1)
r1_spev_abs <- lm(spirit_score_std ~ abs_score_std, d1)

r2_spev_abs <- lm(spev_score_std ~ abs_score_std, d2)
r2_dse_abs <- lm(dse_score_std ~ abs_score_std, d2)

r3_spev_por <- lm(spirit_score_std ~ por_score_std, d3)

r4_spev_por <- lm(spev_score_std ~ por_score_std, d4)
r4_dse_por <- lm(dse_score_std ~ por_score_std, d4)
r4_spev_pv <- lm(spev_score_std ~ pv_score_std, d4)
r4_dse_pv <- lm(dse_score_std ~ pv_score_std, d4)
r4_spev_abs <- lm(spev_score_std ~ abs_score_std, d4)
r4_dse_abs <- lm(dse_score_std ~ abs_score_std, d4)

df_allb %>% 
  full_join(df_all %>% 
              filter(var == "spiritual experience") %>%
              select(-var) %>%
              rename(outcome = scale) %>%
              mutate(predictor = "country")) %>%
  mutate(predictor = factor(predictor,
                            levels = c("country", "porosity vignettes",
                                       "porosity scale", "absorption"))) %>%
  mutate(percent_exp = paste0(round(rsq * 100), "%")) %>%
  select(-rsq) %>%
  spread(study, percent_exp) %>%
  mutate_at(vars(starts_with("study")),
            funs(case_when(is.na(.) ~ ".", 
                           TRUE ~ .))) %>%
  kable() %>% 
  kable_styling() %>%
  collapse_rows(1:2)

Joining, by = c("predictor", "outcome", "study", "rsq")
Column `predictor` joining factor and character vector, coercing into character vectorColumn `outcome` joining factor and character vector, coercing into character vector

predictor	outcome	study 1	study 2	study 3	study 4
country	DSE	.	.	50%	50%
	spiritual events	22%	24%	38%	31%
porosity vignettes	DSE	.	.	.	28%
	spiritual events	13%	.	.	29%
porosity scale	DSE	.	.	.	44%
	spiritual events	.	39%	.	33%
absorption	DSE	.	.	3%	2%
	spiritual events	11%	.	12%	8%

Relationships

All studies, multipart plot

Spiritual Events only, cowplot

fig_s1_por <- d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = pv_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")
# fig_s1_por

fig_s1_abs <- d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s1_abs

fig_s1_title <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s1 <- plot_grid(
  fig_s1_title,
  plot_grid(fig_s1_por, fig_s1_abs, ncol = 1, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))

Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing missing values (geom_point).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing missing values (geom_point).

# fig_s1

fig_s1_title_vert <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s1_vert <- plot_grid(
  fig_s1_title_vert,
  plot_grid(fig_s1_por, fig_s1_abs, ncol = 2, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))

Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing non-finite values (stat_smooth).Removed 20 rows containing missing values (geom_point).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing non-finite values (stat_smooth).Removed 28 rows containing missing values (geom_point).

# fig_s1_vert

fig_s2_abs <- d_all %>%
  filter(study == "Study 2") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s2_abs

fig_s2_title <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s2 <- plot_grid(
  fig_s2_title,
  plot_grid(NULL, fig_s2_abs, ncol = 1, labels = c("", "C")),
  ncol = 1, rel_heights = c(1, 10))

Removed 10 rows containing missing values (geom_smooth).

# fig_s2

fig_s2_title_vert <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s2_vert <- plot_grid(
  fig_s2_title_vert,
  plot_grid(NULL, fig_s2_abs, ncol = 2, labels = c("", "C")),
  ncol = 1, rel_heights = c(1, 10))

Removed 10 rows containing missing values (geom_smooth).

# fig_s2_vert

fig_s3_por <- d_all %>%
  filter(study == "Study 3") %>%
  ggplot(aes(x = por_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")
# fig_s3_por

fig_s3_title <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s3 <- plot_grid(
  fig_s3_title,
  plot_grid(fig_s3_por, NULL, ncol = 1, labels = c("D", "")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s3

fig_s3_title_vert <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s3_vert <- plot_grid(
  fig_s3_title_vert,
  plot_grid(fig_s3_por, NULL, ncol = 2, labels = c("D", "")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s3_vert

fig_s32_vert <- plot_grid(
  plot_grid(
    fig_s3_title_vert, 
    plot_grid(fig_s3_por, labels = c("C")), 
    ncol = 1, rel_heights = c(1, 10)),
  plot_grid(
    fig_s2_title_vert,
    plot_grid(fig_s2_abs, labels = c("D")), 
    ncol = 1, rel_heights = c(1, 10))
)

Removed 10 rows containing missing values (geom_smooth).

fig_s4_por1 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = pv_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_por1

fig_s4_por2 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = por_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_por2

fig_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_abs

fig_s4_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 190))

fig_s4 <- plot_grid(
  fig_s4_title,
  plot_grid(plot_grid(fig_s4_por1, fig_s4_por2, ncol = 2, labels = c("E", "F")), 
            plot_grid(NULL, fig_s4_abs, NULL, ncol = 3, rel_widths = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 1),
  ncol = 1, rel_heights = c(1, 10))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 11 rows containing missing values (geom_smooth).

# fig_s4

fig_s4_title_vert <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s4_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(plot_grid(fig_s4_por1, fig_s4_por2, ncol = 1, labels = c("E", "F")), 
            plot_grid(NULL, fig_s4_abs, NULL, ncol = 1, rel_heights = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 2),
  ncol = 1, rel_heights = c(1, 20))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).Removed 11 rows containing missing values (geom_smooth).

# fig_s4_vert

fig_legend <- get_legend(fig_s4_por1 + theme(legend.position = "bottom"))

Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing non-finite values (stat_smooth).Removed 2 rows containing missing values (geom_point).

fig_all <- plot_grid(fig_s1, fig_s2, fig_s3, fig_s4, ncol = 4,
                     rel_widths = c(1, 1, 1, 2), scale = 0.95)
fig_all

fig_all_vert <- plot_grid(fig_s1_vert, fig_s32_vert, fig_s4_vert, fig_legend,
                          ncol = 1, rel_heights = c(1, 1, 2, 0.2))
fig_all_vert

Daily Spiritual Experiences only, cowplot

fig_s2_abs <- d_all %>%
  filter(study == "Study 2") %>%
  ggplot(aes(x = abs_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s2_abs

fig_s2_title_vert <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s2_vert <- plot_grid(
  fig_s2_title_vert,
  plot_grid(fig_s2_abs, ncol = 1, labels = "B"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s2_vert

fig_s4_por1 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = pv_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_por1

fig_s4_por2 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = por_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_por2

fig_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = abs_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_abs

fig_s4_title_vert <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s4_por1_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_por1, ncol = 1, labels = "A"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_por1_vert

fig_s4_por2_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_por2, ncol = 1, labels = "C"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_por2_vert

fig_s4_abs_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_abs, ncol = 1, labels = "D"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_abs_vert

fig_legend <- get_legend(fig_s4_por1 + theme(legend.position = "bottom"))

fig_all_vert <- plot_grid(
  plot_grid(fig_s4_por1_vert, fig_s2_vert, 
            fig_s4_por2_vert, fig_s4_abs_vert,
            ncol = 2),
  fig_legend,
  ncol = 1, rel_heights = c(2, 0.2))
fig_all_vert

Other versions

Spiritual Events only, one grid, new layout

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  gather(pred_scale, pred_score, c(por_score, pv_score, abs_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption"),
         pred_type = case_when(pred_scale == "Absorption" ~ "Absorption",
                               grepl("Porosity", pred_scale) ~ "Porosity",
                               TRUE ~ NA_character_),
         pred_type = factor(pred_type, levels = c("Porosity", "Absorption")),
         study_scale = paste(study, pred_scale, sep = ": "),
         study_scale2 = case_when(
           study == "Study 4"  & pred_scale != "Absorption" ~ 
             paste(study, pred_scale, sep = ": "),
           TRUE ~ study),
         study_scale3 = case_when(
           study == "Study 4" & pred_scale == "Porosity Scale" ~ "Porosity Scale",
           study %in% c("Study 1", "Study 4") ~ "Porosity Vignettes",
           study == "Study 2" ~ "Porosity Scale", 
           TRUE ~ " "),
         study_scale3 = factor(study_scale3,
                               levels = c("Porosity Vignettes",
                                          "Porosity Scale", " "))) %>%
  filter(!is.na(pred_score),
         spirit_scale == "Spiritual Events") %>%
  ggplot(aes(x = pred_score, y = spirit_score)) +
  facet_grid(rows = vars(pred_type), cols = vars(study, study_scale3)) +
  # facet_grid(pred_type ~ study_scale3) +
  geom_point(data = . %>% distinct(study, study_scale3, country, 
                                   pred_type, pred_scale, pred_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on predictor scale (Porosity Vignettes, Porosity Scale, or Absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Spiritual Events only, by study

fig_s1 <- d_all %>%
  filter(study == "Study 1") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s2 <- d_all %>%
  filter(study == "Study 2") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s3 <- d_all %>%
  filter(study == "Study 3") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s4 <- d_all %>%
  filter(study == "Study 4") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Spiritual Events only, full grid

d_all %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Broken down by predictor type and spiritual scale

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  gather(poros_scale, poros_score, c(por_score, pv_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         poros_scale = recode_factor(poros_scale,
                                     "pv_score" = "Porosity Vignettes",
                                     "por_score" = "Porosity Scale"),
         study_scale = paste(study, poros_scale, sep = ": ")) %>%
  filter(!is.na(poros_score)) %>%
  ggplot(aes(x = poros_score, y = spirit_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                                   poros_scale, poros_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on porosity measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         study_scale = paste(study, "Absorption scale", sep = ": ")) %>%
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = abs_score, y = spirit_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                                   abs_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.2) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on absorption measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")